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ABSTRACT 

An important factor which affects performance of solar adaptive optics (AO) systems 
is the accuracy of tracking an extended object in the wave-front sensor. The accuracy 
of a center of mass approach to image shift measurement depends on the parameters 
applied in thresholding the recorded image, however there exists no analytical predic¬ 
tion for these parameters for extended objects. Motivated by this we present a new 
method for exploring the parameter space of image shift measurement algorithms, 
and apply this to optimise the parameters of the algorithm. Using a thresholded, 
windowed center of mass, we are able to improve centroid accuracy compared to the 
typical parabolic fitting approach by a factor of 3 x in a signal to noise regime typical 
for solar AO. Exploration of the parameters occurs after initial image cross-correlation 
with a reference image, so does not require regeneration of correlation images. The 
results presented employ methods which can be used in real-time to estimate the er¬ 
ror on centroids, allowing the system to use real data to optimise parameters, without 
needing to enter a separate calibration mode. This method can also be applied outside 
of solar AO to any field which requires the tracking of an extended object. 

Key words: instrumentation: adaptive optics, methods: data analysis, techniques: 
image processing, atmospheric effects, sun: granulation 


1 INTRODUCTION 

Tracking extended objects from image sequences in the pres¬ 
ence of noise is required in many different fields. Within as¬ 
tronomy it is used in AO, both for granular images of the sun 
during the day (Scharmer et al.||2003 Rimmele & Radickb 


|1998[ [Michau, Rousset &: Fontanella| 1993|l, a nd elongated 
laser guide stars at night ( Thomas et al.||2008 l. Image shift 
measurement is also used in other fields, for tracking biolog¬ 
ical samples ( [Hand et al.||200^ , and video motion tracking 
applications. 

In solar AO, Shack-Hartmann wave-front sensors 
( Shack &: Platt|[l971 1 with large fields of view are typically 
employed. The cameras used in these sensors have large full 
well depths, as the instruments are photon-noise limited. 
In this work, we investigate tracking extended objects us¬ 
ing data acquired from the Swedish Solar Telescope on-line 
gallery ( Scharmer et al.||1999 l as our wave-front sensor im¬ 
ages (Fig. HF , which have an rms contrast of 10%. We assume 
a photon-noise limited camera with a signal defined as the 
peak intensity above the background, and a noise level de¬ 
fined by the photon-noise. For camera pixels with a typical 
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full well depth of 40000 electrons; the signal would be 4000 
electrons, with photon-noise from 40000 electrons, giving a 
signal to noise ratio (SNR) of 20. 

Shift measurements of extended objects in solar AO are 


calculated in a two-step process (Michau et al. 2006 1 . Ini¬ 


tially an integer shift measurement is performed by locating 
the peak of a cross-correlation of the image with a reference 
image ( Miura et al.|200^ |. Secondly the sub-pixel shift is es¬ 
timated. The determination of the peak location to sub-pixel 
accuracy limits the accuracy to which the shift measurement 
can be performed. 

We concern ourselves with how to best estimate the 
peak location to a sub-pixel accuracy for an arbitrarily 
shaped correlation function derived from cross correlating 
the object with a reference image. For point sources and 
the resultant Airy functions there are analytical methods to 
determine optimal parameters for peak location at a given 
SNR ( Pan, Yang k, Liu||2008 |. However, no such analytical 
treatment exists for images of arbitrary content, such as the 
results of a correlating wave-front sensor. Motivated by this 
we developed a method to optimise the parameters for a 
windowed, thresholded, center of mass measurement for a 
given SNR. We compare this technique with an analytic 2D 
parabolic fit to the central 3x3 pixel region around the cor- 
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Figure 1. Image of solar granulation used as the input image 
in the simulations. The full image is 75 x 75 arc-seconds. Small 
regions of the image are taken and then shifted with respect to 
each-other to artificially generate shifts similar to the effect of 
turbulence of the atmosphere. One such region is shown to the 
right of the full image. It has been re-sampled to the resolution 
used in the simulations, 0.4arc — seconds/pixel. Data obtained 
from the Swedish Solar Telescope On-line Gallery ||Carlsson et al.| 
|2003| |. 


relation peak, as described in Lofdahl (20101. This method 
was chosen as a comparison as its performance is similar 
to the 2D quadratic interpolation method, and significantly 
better than the ID techniques ( L6fdahl|2010 1 and Gaussian 
fitting algorithms ( Waldmann||2007 1. 


2 CORRELATION IMAGE GENERATION 


Simulations for this paper were run on images containing 
solar granulation of a size, field of view and contrast typical 
for solar AO, taken from the large image shown in Fig. 

(Scharmer et al.||1999 1. The image was shifted and binned 


in order to generate images containing sub-pixel shifts using 


the Python language, and numpy routines (Van der wait, 
Colbert &: GaS|2011 1. The solar granulation case used here 


is an example to demonstrate the use of the centroiding tech¬ 
nique, though in general it should be applicable to any ex¬ 
tended image. 

Most of the computational load in centroiding extended 
objects lies in cross-correlating images. By varying parame¬ 
ters applied to centroiding the correlation images, the same 
correlation image can be used. 

Regions of 240 x 240 pixels were taken from Fig.[^ corre¬ 
sponding to a 9.6 arc-second field of view. Integer shifts were 
performed on the full resolution image, with a Gaussian dis¬ 
tribution of mean 0 and standard deviation of 1 pixel, then 
the resultant images were binned by a factor of 10 and had 
shot noise applied, creating typical sub-aperture images of 
24 X 24 pixels, making fully described shifted images down to 
0.1 pixels. These values were chosen to be indicative of resid¬ 
uals in a closed loop AO system. The images, along with the 
known applied shifts were used to compare the windowed 2D 
parabolic fit ( L6fdahl|2010 1 and the windowed, thresholded 
center of mass methods. 


3 PEAK LOCATION ON A CORRELATION 
IMAGE 


3.1 Windowed Parabolic Fit 


A small 3x3 region around the peak of the correlation 
image can be fitted by a 2D parabola, as described in |LofdahT| 
(20101. The parabola takes the form: 


f{x, y) = ai+ a 2 X + a^x^ + a^y + + a^xy, (1) 

where the location of the minima, in x and y respectively 
are given analytically by: 

Xmin — fmin (2^20:5 4:(l2Cl5^ (2) 

ymin — jmin ( 2 ( 23^4 Q,2CLq') /{,^6 4 ^ 305 ^. (3) 


where imin and jmin are the integer positions of the peak of 
the correlation in x and y respectively, and the solution to 
a least squares fit can be found analytically: 

02 = /2 

03 = - 2 {s-i,j)j + /2 

04 = GSi,i). - (Si,-i)j/2 . (4) 

05 = ((Si,l)i - 2 /2 

06 = (si.i - s_l,l - Si,_i -I-/4 


where s describes the 3x3 windowed region around the 
correlation peak, describes the T*' and j**' element of s, 
and i,j can take values from —1 to 1 around the center of 
the peak (located at so,o). 

In high SNR the limiting error in this technique arises 
from the biased sampling of the core of the correlation peak, 
illustrated in Fig. The sampling of the correlation peak 
results in a systematic rounding effect which biases the shift 
estimates towards integer values. The cause of this error is 
apparent in Fig. Here we see the regions windowed for 
use in the centroid highlighted in red. This is a good mask 
for Fig. |3(a)[ however centering on the brightest pixel in 
Fig. 3(b) shows that the peak is being under sampled, and 
not taking into account the full shape of the peak, giving an 
incorrect estimate of the peak location. 


3.2 Windowed, adaptive thresholding Center of 
Mass 

The simplest way to avoid under sampling the correlation 
peak is to use a larger window, however this allows more 
noise into the shift estimate. The noise can be removed to 
some extent by using a threshold to reject contributions from 
parts of the signal, of a similar strength as the noise. For a 
given autocorrelation shape and noise level there will be an 
optimal window size and threshold value, which gives the 
best estimate of the image shift. 

Our proposed method is a two step process. Initially a 
window is placed around the correlation peak, then a thresh¬ 
olded center of mass is taken of the windowed region. The 
size of the window function and the threshold value are vari¬ 
able for each set of images. The threshold value is taken as a 
fraction of the relative peak intensity (max-min of the whole 
correlation image). Re-normalising intensity for every image 
is sympathetic to the shape and size of the correlation peak, 
and ensures that proportionally the same amount of the core 
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Shift Estimation 



Shift applied to image (pixels) 

(a) 



(b) 

Figure 2. Measured image shift plotted against the actual shift 
applied to images. The negative y shifts are plotted in |(a)| to 
make them easier to distinguish. There is a “wobble” apparent in 
the two lines, which is more clearly visible as a systematic effect 
in |(b)[ where the residuals are plotted and take a “sawtooth” 
like pattern. This aliasing effect arises from under-sampling the 
correlation peak. 


of the peak is used in every measurement of the image shift, 
reducing bias effects. 

The correlation image initially has a threshold applied, 
where pixels are rejected if their intensity falls below the 
threshold level, defined by: 

hhresh ^ (fmax imin) ^ PCf, (5) 

where Ithresh is the threshold intensity, I-max is the maximum 
intensity in the correlation image, Imin is the minimum in¬ 
tensity of the correlation image and pet is the fractional 
threshold value. The thresholded correlation image then is 
masked to the chosen window size and is background sub¬ 
tracted, where the background value is the threshold inten¬ 
sity. The centroid estimate of image i, using a reference im¬ 



position (pixels) 

(a) 



position (pixels) 

(b) 

Figure 3. An illustrative ID cut through a correlation peak, 
with the region used by the parabolic fit highlighted in red. Using 
only 3 pixels around the correlation peak, the shift estimate can 
be unavoidably biased away from the true location of the peak. 
While [(a)] shows the ideal case for using this method, there are 
some cases where the shift differs from the measured position due 
to the limited size of the region used, as demonstrated in |(b)[ This 
is shown by the arrows above the plots, the green arrow indicates 
where the parabolic centroid estimates the correlation peak, while 
the blue arrow shows the true location of the peak. 


age r, can be described as a vector 


R 


i,r 


*0 

Vo 


( 6 ) 


where xq and yo are the x and y components of the centroid 
estimate R*’”. R*’” is calculated using: 

2/max fCmax 

1 ^ ^ (7) 

y=l x=l 


where I is the total intensity of the correlation image, Ix,y is 
the intensity of pixel x,y in the correlation image with the 
threshold applied, and R^’.y is the vector position of [x, y] in 
the correlation image. 

The size of the window is a relatively small parame¬ 
ter space to explore, going from a single pixel around the 
core (corresponding to an integer shift measurement), to 
the wings of the correlation peak succumbing to background 
noise. If any larger boxes are used a drop off in performance 
is seen as more noise is included in the centroid estimate, 
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without any extra useful information being added. The outer 
threshold for this parameter needs to be set arbitrarily. If 
too small a window is used, a similar effect to the windowed 
parabolic fit is seen, in that the measurements are biased 
toward integer shifts. The optimum window size is chosen 
as a trade off between including as much of the correlation 
peak as possible, but also minimizing the number of pixels 
which only contribute noise to the measurement. 

The centroid threshold value is normalised such that a 
value of 1 uses only pixels with the maximum flux and a 
threshold of 0 uses all available pixels. This parameter be¬ 
haves similarly to the window size, in that using more pixels 
increases the noise contribution, reducing the accuracy of the 
shift estimate. Using high thresholds gives rise to a bias to¬ 
wards integer shift measurements, similar to that seen in the 
parabolic fit. The optimum threshold value lies somewhere 
between these two regimes, and is liable to change depend¬ 
ing on the window size. This means the whole parameter 
space needs to be explored for all window sizes to identify 
the best combination of parameters for the centroids. We 
optimise the threshold and window size for a set of images, 
which all use a common reference image. The parameters 
then only need to be updated when the reference image is 
changed to take into account slow changing effects, such as 
the evolving granulation pattern on the solar surface. 

The optimum set of parameters will depend on a num¬ 
ber of different obvious factors, including the image shape, 
the shape of the resulting correlation function, and the SNR 
of the images, assuming an arbitrary unknown correlation 
shape. There is no obvious analytic way to determine the 
best parameters for a given set of images, or circumstances, 
hence we explore the parameter space to find the optimal 
solution. However once the optimum set of parameters is 
found for a given object, at a set SNR level, then it should 
be constant until one of these factors changes. In solar AO 
the regions used for wavefront sensing are constantly evolv¬ 
ing, causing the reference image used to be updated on a 
frame by frame basis. This also means that over time the 
optimum parameters are subject to change and need to be 
updated. As the parameters chosen are based on normalised 
intensity, they are insensitive to changes in flux for a given 
SNR, such as scintillation effects. 


3.3 Error estimation 

Given a set of shifted images of matching content and SNR, 
it is possible to make multiple independent estimates of the 
image shift. By comparing the spread of the shift estimates 
we can get an estimate of the error on the shift measure¬ 
ment. Using different reference images allows us to estimate 
the shift in the image multiple times, we can then use the 
standard deviation of the shift estimates as an indicator of 
the error on the shift estimate. As this is a statistical pro¬ 
cess, the estimated error will not be accurate for the shift 
estimate of a single image in the set. However when aver¬ 
aged over the set of images, we can estimate the magnitude 
of the shift error on the set. 

This set of images may be drawn from a single temporal 
wave-front sensor frame in AO, guaranteeing spatial similar¬ 
ity of the images. Alternatively the set could be drawn from 
a time sequence in correlation video tracking. Care must be 


taken that the object does not change its spatial character¬ 
istics significantly over the duration of the set. 

We use multiple different reference images, e.g. using 
the first 10 sub-aperture images in the wave-front sensor 
frame, and use each of them as a reference to estimate image 
shifts. The global tip/tilt terms are then removed to compen¬ 
sate for the systematic error in shift estimation, due to the 
unknown shift applied to the reference image. This is a com¬ 
mon practice in AO systems to negate effects like wind shake 
from measurements. The subtraction of the global tip/tilt 
term can be described with: 

= R’- - (R'-), , (8) 

for a given reference image, where R^/j is the center of mass 
estimate of a set of images with tip/tilt removed, and (R"^),, 
is the averaged tip/tilt term over all of the images using a 
given reference. This removes the shift due to each of the ref¬ 
erence images, making the centroid estimates from different 
reference images directly comparable. The standard devia¬ 
tion of the resultant shifts estimate the error, . This 

t/t 

method of estimating centroiding errors allows for the pa¬ 
rameter space to be explored on real data, where the actual 
shifts are unknown, and not just on simulated data. 


4 RESULTS 


The full parameter space was explored in simulation for a 
range of threshold values and window sizes applied to the 
correlation images. Fig. |4(a)] shows the magnitude of residual 
errors for different sets of parameters. Fig. |4(b) shows the 
standard deviation of the centroid measurements using ten 
different reference images. This has the same characteristics 
as the real error values, showing it can be used to estimate 
the location of the optimum parameters for the centroiding 
algorithm. The optimal parameters from each of the meth¬ 
ods is highlighted with a white marker. 

In the thresholding axis (x) of Fig. I^it is possible to see 
the effects of aliasing towards the large thresholding values 
on the right of the plots. This effect is similar to the aliasing 
in the parabolic fit, and in all cases the error approaches 
that of integer pixel estimation, as at the largest threshold 
value only the brightest pixel is considered, equivalent to an 
integer pixel shift estimation. 

In the window size axis (y) of Fig. the structure is 
more complicated. Initially the aliasing is apparent for small 
window sizes, similar to the parabolic fit. This problem de¬ 
creases as the window size increases, until its optimal region. 
However the performance begins to degrade again for large 
windows for low thresholding values. This happens where 
the region is so large that as well as including all of the peak 
of the correlation, it includes increasing amount of noise, 
which isn’t filtered out by the thresholding. 

The centroid optimisation was performed on a range 
of different noise levels (using photon-noise) to demonstrate 
how noise affects the centroid estimates. The parameters de¬ 
pendence on SNR is demonstrated in Fig. with Fig. 5(a) 
showing how the threshold level affects the accuracy of the 


centroid estimates, and Fig. 5(b) illustrating how changing 
the window size affects the accuracy of the centroid esti¬ 
mates. The estimation was performed on ten different re- 
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Figure 4. Full parameter space for the box size and threshold 
value in the center of mass algorithm. |(a)| shows the real error 
associated with the parameters used in the center of mass tech¬ 
nique, and 1(b) I gives the error estimate taken from the standard 
deviation on centroids using multiple reference images. The shape 
of the two plots is similar, indicating multiple references is a suit¬ 
able estimator of the error. The white spots on the plot show 
where the optimum parameters lie for the respective methods. 
The estimated error position does not directly overlap with the 
location of the real minima, but it can be seen that the difference 
in error is minimal. 


gions of the granule image, with the errors taken to be the 
standard error. 

The optimal values for the parameters varied with SNR 
as can be seen in Fig. Fig. 6(a) shows the optimal thresh¬ 
olding values for the various SNR levels, both best perform¬ 
ing and the best estimated. The estimated threshold levels 
differ from the true value up to a SNR of 2, where the esti- 



(a) 



Figure 5. |(a)| shows how the optimum threshold value is affected 
by different SNR levels, [(^ shows how the window size affects the 
error on the centroid estimate for different SNR levels. Above a 
SNR of 5 the curves no longer change, staying at their high SNR 
shapes. 


mated threshold value is consistent, and in a region where 
small variations have little effect on the accuracy of the shift 
estimate. This trend is also seen in Fig. 6(b)[ at low SNR lev¬ 
els the estimated box size is larger than the actual optimal 
value, but at higher SNR levels they agree more. 

The optimal parameters for thresholding and window 
generally reduce the number of pixels used in the centroid 
in low SNR conditions by using small thresholds and small 
window sizes to reduce the amount of noise in the centroid. 
At higher SNR values the parameters stabilise for a given 
set of images to give the most accurate shift estimate. 

Our technique fails in the low SNR regime. This is due 
to the different sources of noise in the correlation image, and 
our sampling of them. The simplest way to do this is to have 
the images and their noise terms separate, as in equation]^ 


Im — ImsjgTidZ T ImTtoise 

Ref = Ref signal Ref noise 


( 9 ) 


where Im represents the overall image being centroided, 
Imsignai describes the signal in the image, and Imnoise de¬ 
scribes the noise associated with the image, in our case shot 
noise. Ref follows similar definitions for the reference image. 
When combined, assuming a linear regime, the correlation 
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(b) 


Figure 6. |(a)| shows the optimal thresholding value for the dif- 
ferent SNRs of the images used in the centroiding. Initially the 
thresholding is high, to remove as much noise as possible from 
the correlation image, then the thresholding drops to its optimum 
value for images which have low noise. [(^ shows the box size for 
the different SNR levels. This shows a similar tend, of increas¬ 
ing window size at high SNR, using more pixels when the noise 
is reduced. At low SNR the estimated parameters disagree with 
the true optimal parameters, but this disagreement decreases at 
higher SNR. 


image has four terms: 

Corr = Corr + Corr 

Corr [^Im^oiseRefs^g^g;] “t" Corr [imjio^seRef^oise] (10) 

where Corr is the total signal in the correlation image, with 
the contributing factors all described to the right. If we as¬ 
sume that the contribution of Corr [im„oi 3 gRef„oige] is negli¬ 
gible, then there are two remaining error terms which affect 
our estimate of the centroid. However by taking an average 
over different references in our estimate of the error, we are 
in effect averaging out the Corr [imgig„g,Ref„gije] term. This 
term becomes more dominant at lower SNR levels, hindering 
the performance of our technique. There are other methods 
of estimating the error of a centroid on an extended object. 


such as Saunter (20101, which don’t have this problem, but 



Figure 7. This plot shows the performance of the center of mass 
algorithms and the 2D parabolic fit for a range of different SNRs. 
It can be seen above a SNR of 1 the windowed, thresholded cen¬ 
ter of mass outperfroms a 2D parabolic fit. The 2D parabolic fit 
tapers off in performance at 0.05 pixel error, whereas the win¬ 
dowed center of mass has a much lower performance threshold. 
The vertical line on the plot show the expected SNR for a solar 
granule image with a contrast of 10%, and a camera with a full 
well depth of 40000 electrons, which represents typical conditions 
in solar AO. It can also be seen that the performance from esti¬ 
mating the errors on the center of mass is worse than the optimal 
case, but does still reach close to peak performance. 


this requires an oversampling of the correlation peak, some¬ 
thing avoided in AO to reduce data rates and computation 
time. 

The overall performance of the centoriding techniques 
for the different SNRs is shown in Fig. For high SNR, 
the best performance is given by the thresholded, windowed 
center of mass measurements, with little difference between 
the theoretical best performance and the performance de¬ 
rived from error estimation. The overall boost in accuracy 
is 3x for the high SNR. For SNR below 1, the windowed 
parabolic fit outperforms the thresholded, windowed center 
of mass method. This could be due to the crude error esti¬ 
mator implemented here and it may be possible to improve 
this using other error estimation techniques ( |Saunter|2010p . 
However this still could only bring the performance back to 
the level of the 2D parabolic fit at best. Our technique is 
best suited to high SNR regimes. 


5 CONCLUSIONS 

We have demonstrated that for tracking extended sources, a 
method of error estimation allows different centroiding pa¬ 
rameters to be explored on real data, allowing for the opti¬ 
mum parameters to be chosen. While this does take extra 
computation, the correlation images only needs to be gen¬ 
erated once for each reference, minimizing the increase in 
computation effort required. Also once the optimum set of 
parameters has been found, they should hold as the best pa¬ 
rameters until something in the system changes, i.e. a change 
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of target, or reference image. The parameters for the cen- 
troiding algorithm should be updated regularly to keep it 
optimal. 

Exploring the parameter space is a parallelisable pro¬ 
cess, so can be performed quickly. With the use of SIMD 
(Furht 20081 and more advanced optimisation algorithms, 
rather than the brute force method exploring the full param¬ 
eter space implemented here, the method should be viable 
for use in a real-time system. 

The method of noise estimation used here is crude, 
though good enough for our purposes, and could be used for 
different parameters in other techniques, such as |Li et alT] 
(|2008|). There are more efficient algorithms for estimating 
noise on centroiding of extended objects, such as |Saunter| 
(20101, which could also be implemented to give quanti¬ 
tative estimators of centroiding accuracy, as well as being 
computationally less intensive than the multiple reference 
approach. 

Overall, for the solar case, with high SNR, the use of an 
optimised, thresholded, windowed center of mass algorithm 
offers a factor of 3x improvement in centroiding accuracy 
over the windowed parabolic fit. This could be used real¬ 
time in solar AO for better wave-front estimations, and also 
with post processing techniques, such as measuring more 
accurate atmospheric profiles. 

Further investigation should be performed in the low 
SNR regime, where both the center of mass and 2D parabolic 
fit methods give poor performance, to see if more accurate 
centroids can be extracted. There is also more work to be 
done in implementing the technique into a real system which 
performs centroiding on extended objects, to see how it af¬ 
fects system performance. 
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